Fast Pixel Space Convolution for CMB Surveys with 
Asymmetric Beams and Complex Scan Strategies: FEBeCoP 

S. Mitra 
Sanj it . MitraOjpl . nasa . gov 
^^ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA 

o 

(N 

-J G. Rocha 

^ graca@caltech.edu 

^^ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA 

<N 
<N 

I I K. M. Gorski^ 

^^ Krzysztof . M . Gorski® jpl . nasa . gov 

\^ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA 

• 

I K. M. HuflFenbercer^ 

O 

Y^ huf f enbeOphysics.micLmi .edu 

^ Department of Physics, University of Miami, 1320 Campo Sano Drive, Coral Gables, FL 33146, USA 

^ H. K. Eriksen^ 

0^ h.k.k.eriksen@astro.uio.no 



(N 



^ 



Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway 



g M. A. J. Ashdown^ 

^D majalOmrao. cam. ac.uk 



Astrophysics Group, Cavendish Laboratory, J. J. Thomson Avenue, Cambridge CBS OHE, UK 

C. R. Lawrence 
Charles .R.Lawrence@jpl.nasa.gov 
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA 



-•^Warsaw University Observatory, Aleje Ujazdowskie 4, 00-478 Warszawa, Poland 

^Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA 

^Centre of Mathematics for Applications, University of Oslo, P.O. Box 1053, Blindern, N-0316 Oslo, Norway 

^Institute of Astronomy, Madingley Road, Cambridge CBS OHA, UK 



ABSTRACT 

Precise measurement of the angular power spectrum of the Cosmic Microwave Background 
(CMB) temperature and polarization anisotropy can tightly constrain many cosmological models 
and parameters. However, accurate measurements can only be realized in practice provided all 
major systematic effects have been taken into account. Beam asymmetry, coupled with the scan 
strategy, is a major source of systematic error in scanning CMB experiments such as Planck, the 
focus of our current interest. We envision Monte Carlo methods to rigorously study and account 
for the systematic effect of beams in CMB analysis. Toward that goal, we have developed a fast 
pixel space convolution method that can simulate sky maps observed by a scanning instrument, 
taking into account real beam shapes and scan strategy. The essence is to pre-compute the 
"effective beams" using a computer code, "Fast Effective Beam Convolution in Pixel space" 
(FEBeCoP), that we have developed for the Planck mission. The code computes effective beams 
given the focal plane beam characteristics of the Planck instrument and the full history of actual 
satellite pointing, and performs very fast convolution of sky signals using the effective beams. In 
this paper, we describe the algorithm and the computational scheme that has been implemented. 
We also outline a few applications of the effective beams in the precision analysis of Planck 
data, for characterizing the CMB anisotropy and for detecting and measuring properties of point 
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1. Introduction 

Cosmic Microwave Background (CMB) anisotropy 
measurements have played a leading role in preci- 
sion cosmology. Detection of the CMB monopole 
by Penzias and Wilson (Penzias & Wilson 1965) 
firmly established that the Big Bang model is 
consistent with observation. The Cosmic Back- 
ground Explorer (COBE) (Smoot et al. 1992) de- 
tected the tiny ~ 10~^ fluctuations, which justi- 
fied how structures could form in a statistically 
isotropic universe starting from seed density per- 
turbations. The Wilkinson Microwave Anisotropy 
Probe (WMAP) (Larson et al. 2010) measured 
cosmological parameters — e.g., the flatness of the 
universe — with unprecedented accuracy. There 
have been many suborbital experiments (E. Tor- 
bet et al 1999; A. D. Miller et al 1999; S. Padin 
et al 2001; T. J. Pearson et al 2003; C. L. Kuo 
et al 2004; P. de Bernardis et al 2000; S. Hanany 
et al 2000) that have probed specific aspects of 
CMB anisotropy. Yet only a small fraction of the 
information embedded in the CMB temperature 
anisotropy has so far been measured, and mea- 
surements of the CMB polarization anisotropy 
have just begun. Precise measurement of the 
CMB temperature and polarization anisotropy is 
the most important next step in building a rig- 
orous model of the content and evolution of the 
universe. The Planck mission is a substantial 



step forward in CMB anisotropy measurement — it 
aims to measure almost all the information avail- 
able in temperature anisotropy and will make a 
full-sky, high-resolution, low-noise measurement 
of CMB polarization anisotropy. It will also im- 
age thousands of point sources of astrophysical 
and cosmological interest. 

Planck (The Planck Collaboration 2006; Tauber 
et al. 2009) is a full-sky scanning experiment 
launched in May 2009, taking data since August 
2009. It will produce high-resolution maps with 
tens of millions of pixels. The Low Frequency 
Instrument (LFI) covers 30, 44, and 70 GHz; the 
High Frequency Instrument (HFI) covers 100, 143, 
216, 353, 545, and 857 GHz. From the second 
Lagrangian point of the Earth-Sun system (I/2), 
Planck scans nearly great circles on the sky, cov- 
ering the full sky twice in one year ((jL^npnC & 
Tauber 2005)). The satellite spins at 1 rpm around 
an axis that is re-pointed roughly 30 times per day 
along a cycloidal path, with the spin axis moving 
in a 7? 5 circle around the anti-Sun direction with 
a period of six months. This ensures that all feeds 
completely cover the sky, including the ecliptic 
polar regions. 

High-resolution, low-noise measurement of the 
CMB anisotropy is challenging. Cosmic-variance- 
limited measurement at high resolution demands 



not only precision, but very high accuracy — ah 
the systematic effects must be accounted for; oth- 
erwise systematic errors can become significant 
or dominant compared to statistical errors. The 
effect of asymmetric beams is one of the most 
important sources of systematic error in current 
and future CMB missions, which invest a signifi- 
cant effort in measuring the polarization. Prob- 
lems with asymmetric beam systematics, which 
could transfer power from large scales to small, 
could be acute for measurements of the rapidly 
falling power spectra at high resolution. For po- 
larization (4shdown et al. 2009), where the signal 
is an order of magnitude below the temperature 
anisotropy and the inherent directionality of the 
polarization field gets coupled to the asymmetry 
of the beam, the effect could lead to large devi- 
ations in the polarization anisotropy. In general, 
precision measurement of the polarization power 
spectrum at high-resolution requires a correction 
for asymmetric beams. 

Planck beams range from ~ 30' at 30 GHz to 
~ 5' at 217GHz and higher frequencies. The low- 
est frequency channels are farthest from the cen- 
ter of the focal plane and exhibit the most asym- 
metric beams. For example, the 30 GHz chan- 
nel has beams with ellipticities of the order 1.37. 
In Figure 1 we show realistic Planck beams, for 
one 30 GHz and one 143 GHz detector, simulated 
by a full diffraction analysis of the telescope us- 
ing GRASP9 (Sandri et al. 2002a,b; Maffei et al. 
2010; Yurchenko et al. 2004), while in Figure 2 
we show the focal plane layout (Dupac & Tauber 
2005). The related problems of measurement of 
beam shapes (nunenoerger et ai. zuiu) and the 
incorporation of their uncertainties in the anal- 
ysis (Rocha et al. 2010) are also active areas of 
research. Our focus in this paper is to study the 
effect of asymmetries in the beam shapes on the 
sky maps. 

Broadly speaking, the problem of asymmetric 
beams at high resolution has been addressed in 
the literature in two ways: estimating the bias 
of the "pseudo-C^" estimator semi- analytically in 
the spherical harmonic basis (Mitra et al. 2009; 
Hinshaw et al. 2007; Mitra et al. 2004; Fosalba 
et al. 2002; Souradeep & Ratra 2001); or decon- 
volution of observed maps (Armit age- C apian & 
Wandelt 2009; Burigana & Saez 2003). However, 
to handle a complicated instrument like Planck, 
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Fig. 1. — Logarithmic plots of Planck beams for 
30 GHz (top) and 143 GHz (bottom) as simulated 
by GRASP9. The 30 GHz beam was simulated in 
Cartesian coordinate and the 143 GHz one in polar 
coordinates. (The difference in coordinate system 
is just a matter of representation, it is not related 
to the detectors or beams.) 



with non-trivial scan strategy and noise character- 
istics, we aim to use a simple forward approach, 
i.e., to estimate the systematic effect by study- 
ing Monte- Carlo- simulated maps convolved with 
detector beams using the real scan path. This re- 
quires a fast convolution algorithm for scanning 
CMB experiments like Planck. Existing simu- 
lation methods (Reinecke et al. 2006), involving 
spherical harmonic based convolution algorithms, 
are relatively slow, and even computationally pro- 
hibitive for a large number of simulations (typi- 
cally a few thousand in Monte Carlo studies). This 
is because they generate time ordered data (TOD) 
prior to mapmaking. Improving simulation capa- 
bilities was the main motivation for developing a 
fast, easy-to-use map convolution method which 
would include the effect of asymmetric beams and 
scan. 

Convolution in pixel space is fast and flexible 
because in scanning CMB experiments the main 
beam is highly localized in pixel space. The stan- 
dard choice of pixel number A^pix for a given beam 
size ensures that only A/'beam ^ 200 pixels are re- 
quired to accurately describe the beam. This is 
much smaller than the number of spherical har- 
monic moments needed to describe the beam with 
adequate accuracy. The total computation cost 
for pixel space convolution (after effective beam 
computation) scales as A/'pix A/'beam • This is very 
small compared to any other existing method. For 
example, using convolution in the harmonic basis 
for an idealized scan (where each direction has a 
fixed beam orientation angle), performed by the 
fast multipole based method "total convolution" 
(Wandelt & vjorsKi zuui; ^^namnor ei ai. zuuu), 
scales as ^ ^max^^max, where ^max is the highest 
multipole moment that contains non-negligible in- 
formation about the anisotropics in the map and 
^max is the highest azimuthal multipole needed 
to adequately describe beam asymmetry. This 
is more demanding than the pixel space convo- 
lution. There has been recent progress in making 
this step faster using recursive algebra (Prezeau & 
ReinecKe zu ). In practice, though, convolution 
is much more complicated than the ideal case, as 
the scan can take the beam to any arbitrary orien- 
tation and one has to go through the whole scan 
path and interpolate between orientation angles 
obtained from the fixed Total Convolution grid, 
which is a tedious process. Considering all these is- 



sues, pixel space is often a more convenient choice. 

We have developed a HEALPix^ (Gorski et ai. 
2005) based pixel space convolution method by 
pre-computing the "effective beams" for each pixel 
for both temperature and polarization maps, in- 
corporating asymmetric beams of the instruments 
and the real scan strategy. The effective beam 
for a given pixel is essentially the average of (dis- 
cretized) beam functions for every pointing sample 
that fall in that pixel (a formal mathematical defi- 
nition is provided later in the paper). The method 
has been successfully implemented for the Planck 
mission as a computer code FEBeCoP, acronym of 
"Fast Effective Beam Convolution in Pixel space" . 
FEBeCoP efficiently utilizes large distributed com- 
puting facilities, hence it is fast, but provides a 
very easy interface, which enables external appli- 
cations to use this facility in a convenient way. 
Different analyses, which require realistic Monte 
Carlo simulations of observed sky maps, but not 
TODs, can benefit from this method. 

In scanning experiments which target static ob- 
jects (not transient sources), the data products are 
better understood by studying the effective beam 
rather than the detector intrinsic beam function. 
For instance, consider an elliptic beam function 
that moves on the sky in a scanning pattern simi- 
lar to that of Planck. If the scan is along the mi- 
nor axis of the ellipse, repeated hits in a pixel will 
tend to make the effective beam more symmetric 
than the true beam, hence the final map will have 
less severe high resolution systematics. On the 
other hand, if the scan is along the major axis, 
the elongation will be enhanced in the effective 
beam and the final map will have more high reso- 
lution systematics. Effective beams thus provide a 
more intuitive picture of the systematic effects in 
the observed map as compared to the detector in- 
trinsic beam function. This feature would be use- 
ful not only to understand the non-trivial effects 
of beams on observed CMB anisotropy, but also 
for precise measurement and extraction of point 
sources and Sunyaev-Zel'dovich (SZ) clusters from 
Planck data (ihe Planck Collaboration 2UUt)) us- 
ing the accurate knowledge of the point spread 
function (PSF) obtained from the effective beams. 

The paper is organized as follows. In sec- 
tion 2 we introduce effective beams. Correspond- 
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ing mathematics for temperature-only and polar- 
ized effective beams are given in section 3. Sec- 
tion 4 analyzes the accuracy of the method com- 
paring with theoretically exact results. Section 5 
describes the details of computational implemen- 
tation, which is followed by two types of applica- 
tions of the effective beams. Simulation of maps 
for Planck and comparison with existing simula- 
tions are shown in section 6. Study of collective 
effect of the beams via elliptical Gaussian fits and 
estimation of approximate scalar transfer function 
for pseudo-C^ estimator by direct spherical har- 
monic decomposition of beams are covered in sec- 
tion 7. We conclude with summary and discus- 
sions in section 8. 

2. Effective Beam: Method & Advantages 

The finite angular response of a detector is de- 
scribed by the beam function, which is closely re- 
lated to the point spread function (PSF) in im- 
age space. The effective beam is the average of all 
beam functions pointing at a certain direction for 
a given scan strategy. Since, in practice, a scan 
strategy may not pass through exactly the same 
point twice, the effective beam concept makes 
sense only when the sky is discretized, then one 
can define an effective beam for each pixel on the 
sky. 

The mapmaking algorithms (Ashdown et al. 
2007; Hinshaw et al. 2007) commonly used in CMB 
analyses do not take into account the finite size of 
the instrumental beams, so the resulting map car- 
ries the signature of the instrumental beam func- 
tion and the satellite scan strategy. In fact, the 
observed map is a convolution of the true CMB 
sky with the effective beams. Effective beams can 
therefore be treated as the primary object for 
studying the combined systematic effects of beams 
and scan strategy. To appreciate this fact con- 
sider the hypothetical case of symmetric detector 
intrinsic beam. Since the beam center (pointing 
direction) need not coincide with the pixel center, 
if the distribution of pointing directions in a pixel 
is elongated (say only one set of scan rings passes 
though that pixel), effective beams can become 
asymmetric. 

If the effective beams were symmetric, high res- 
olution CMB data analysis would be computation- 
ally trivial. However, in practice such an advan- 



tage is absent. The effective beams in CMB exper- 
iments are (significantly) asymmetric due to sev- 
eral unavoidable reasons, which include: 

• intrinsic optics of the instruments 

• not all detectors can be placed close to the 
principal optic axis, so there are aberrations 

• a sample is an integral along the scan path 

• pointing directions are not uniformly dis- 
tributed in a pixel. 

Moreover, the shape of effective beams varies 
across the sky. Our aim is to estimate the ef- 
fect of asymmetric effective beams and account 
for it in data analysis. 

One of the immediate advantage of effective 
beams is that the point spread functions, the 
transpose of effective beams, become easily deriv- 
able at every pixels which is otherwise non-trivial. 
For instance, in order to compute PSF at neigh- 
boring pixels using point source convolution meth- 
ods, convolutions must be performed separately 
for each pixel to prevent overlap between PSFs. 
Limitation of the effective beams method is that 
the PSFs are centered at the pixel centers, though 
it may not be an issue in (pixelized) map based 
applications. 

As mentioned before, Monte Carlo simulations 
are extremely important for studying systematic 
effects. Consider a sky map observed with certain 
detector response and scanning. To study the sys- 
tematics in this map, one would simulate a number 
of model input maps and convolve them with the 
same detector response and scanning to estimate 
the average effect on these simulated maps. This 
can then be used to draw inferences on the map ob- 
served by the instrument. Since at each simulation 
scanning and beam are the same, huge amount of 
computation can be saved by pre-computing the 
effective beams, which include the effect of beam 
and scan once and for all for the whole set of sim- 
ulations. 

Once the effective beam is computed, it is 
rather straightforward and notably fast to make 
simulated "observed" maps. This is indeed the 
main goal of this paper - to show that the effec- 
tive beam convolution is fast and accurate. In the 
end, we also provide some insights on the observed 



systematic effects directly by studying the proper- 
ties of the effective beams. 

The cost of effective beam computation scales 
as A/tod X ^beam, whcrc A/tod is the total number 
of samples, and A^beam is the number of pixels re- 
quired to adequately describe the beam, and cost 
of convolution scales as A^pix x TVbeam , where A^pix 
is the total number of equal area pixels in the map. 
In general, A/'beam = A^pix, as the far side lobe of 
the detector response can span the whole sky. In 
practice, however, the beam is negligible at most 
of the pixels; in fact, for CMB experiments of in- 
terest, non- negligible part of a main beam is very 
small and localized^, leading to a big saving in 
computation and Input-Output (I/O) handling. 

The finite size of the instrumental beam re- 
sponse causes exponential damping of angular 
power spectra at high multipoles, where noise 
becomes dominant. As a result, a CMB map 
observed with a certain Full Width at Half Maxi- 
mum (FWHM), FWHMbeam, has usable informa- 
tion content typically up to a multipole of ^max ^ 
4/FWHMbeam. Also, a CMB map is pixelated in 
such a way that the number of pixels is minimized 
while retaining all non-negligible information. Ac- 
cording to the HEALPix guidelines, this can be 
achieved by choosing A^pix - 200/FWHMbeam- 
HEALPix divides the sky into A^pix = 12 n^^^^ 
pixels, where ngide is an even power of 2 (Gorski 
et al. 2005). Thus the above choice of pixelization 
translates to, ngide ~ -^max = 4/FWHMbeam. As 
an example, a map observed with ~ 30' beam is 
usually tessehated with N^i^ = 3145728 HEALPix 
pixels (i.e., ngide = 512). Then, if we assume that 
the beam is significant only up to an axial distance 
of r times FWHM, the number of pixels required 
to describe the beam would be 



A^be 



47r 
200 



47rFWHM 



50 r^ 



X area (2-1) 

^(r FWHMbeam)' 

(2-2) 



beam 



which means that, if the beam is significant up to a 
diameter of four times FWHM, the number of pix- 
els needed to describe it with adequate accuracy 



^Here we treat the far side lobes separately — they can be 
better modeled in the spherical harmonic domain at a much 
lower angular resolution (Wandelt &; Gorski 2001). 



is ~ 200. Note that, a symmetric Gaussian beam 
falls to 2-16 ^ 0.0015% (-48dB) of its peak value 
at this diameter, and, in section 4.1, we explicitly 
show that the net power (effective area) outside 
this diameter is ~ 0.001% of the total beam power. 

We compare the computational efficiency of 
pixel space convolution with the fast harmonic 
space convolution algorithm "Total Convolution" 
(Wandelt & Gorski 2001; Chahinor et al. 2000). 
Total Convolution uses a Fast Fourier Transform 
(FFT) on the sphere and convolution in the conju- 
gate space to generate a "Total Convolution Data 
Cube" (TCDC) — sets of convolved sky maps for 
an equi-spaced grid of position and orientation co- 
ordinates of the beam. Total Convolution is the- 
oretically exact at those grid points for band lim- 
ited beam functions. Even though Total Convolu- 
tion itself does not incorporate satellite scan strat- 
egy, the algorithm sets the limit on computation 
cost and accuracy that can be achieved by har- 
monic space simulation methods. Below we pro- 
vide approximate theoretically estimated compu- 
tation cost scaling of effective beam and total con- 
volution; later we provide numerical comparison 
of accuracy. Note, however, that effective beam 
method also includes scan strategy (in the pre- 
computation stage), to do the same in harmonic 
space one goes through the whole scan strategy 
and interpolates a TCDC, adding a computation 
cost which scales poorly due to machine I/O. 

Computation cost of Total Convolution scales 
^s "^ ^max^max (Waudclt & Gorski 2001; Challi- 
nor et al. 2000; Prezeau ^ Reinecke 2010), 
where ^max is the highest multipole moment that 
contains non-negligible information about the 
anisotropics in the map and mmax is the high- 
est azimuthal multipole needed to adequately de- 
scribe beam asymmetry. In order to maintain 
minimum accuracy in convolution, one needs to 
use Cax ^ \/^^pix/3 (= 2nside) and m^ax > 10, 
which means that the computational cost would 
scale as > nJ^ />/3. Since convolution in pixel 
space scales as ~ A^pix^beam, pixel based convo- 
lution would be faster than total convolution as 
long as A^pix > 3 x N^^^^ -- 10^, as A/'beam was 
shown to be few hundred for a reasonable beam 
cutoff radius. 

As effective beam computation and convolution 
are the focus of the paper, the rest of the pa- 
per provides details of computation and accuracy 



analysis considering the main beam only. In addi- 
tion, though our method and implementation can 
handle any arbitrary beam shape, in order to reli- 
ably compare our results with existing Planck sim- 
ulations we consider only elliptic Gaussian beam 
for numerical results. 

3. Mathematics for Effective Beams 

Computation of effective beams involves large 
number of mathematical operations. Hence, it is 
important to formulate it in a computationally 
minimalist manner. In this section we list the for- 
mulae we use in our code. Before we proceed fur- 
ther we define some quantities which will be used 
throughout the paper: 

• Effective beam for pointing pixel i and beam 
pixel j: Bij = B. 

• Pointing matrix for pixel p and sample (or 
time) t: 



A 



tp 



1 , if sample lies in pixel 
, otherwise 

(3-1) 



• Measured signal or Time Ordered Data 
(TOD) at sample t: dt = d. 

• Pointing direction at a sample t: f(t) = 

• Polarization angle, the angle between the po- 
larization axis of the detector and the local 
meridian: t/jt = ip{t). 

(for unpolarized detectors, polarization axis 
is the X-axis of the fiducial detector frame) 

• The set of three pointing angles that define 
the instantaneous geometry of each detector 

• Angle between polarization axis and the x- 
axis of the fiducial detector frame: t/jpoi. 

• Angle between the major axis of an elliptical 
beam (or any preferred axis relative to which 
the beam profile has been defined) and the 
X-axis of the fiducial detector frame: '^eii- 

• Detector intrinsic beam function at a direc- 
tion 1*5 = [^5, ^5] with the pointing angles p: 

^(r6,p). 



• Direction vector corresponding to pixel num- 
ber p: Up. 

Note that the polarization angle is being used 
to describe the orientation of the detector with- 
out any loss of generality. Since the polarization 
axis of the detector is fixed to the fiducial detector 
frame, any other orientation angle in this frame, 
e.g., the orientation angle of the major axis of an 
elliptic beam, '^eii, can be derived by adjusting the 
appropriate offset angles. For an unpolarized de- 
tector, the angle between the polarization axis and 
the fiducial detector frame is taken as zero. 

The detector intrinsic beam function is gener- 
ally given in a 2-D coordinate system on the beam 
plane, with the origin at the beam center and 
the X-axis conveniently chosen — e.g., for an ellip- 
tic Gaussian beam the major axis is chosen as the 
X-axis. Here we consider the beam to be defined 
on the tangent plane, in general the beam could as 
well be specified in the "curved" polar coordinates 
at the North pole, then we could take the = 
direction as the x-axis. 

3.1. Temperature- Only Effective Beam 

By definition, effective beams are the objects 
whose convolution with the true CMB sky should 
produce the observed sky map. So we write the 
expression for the observed CMB map (without 
noise) in a convolution form and then define effec- 
tive beams. 

Each observed sample is a convolution of beam 
and the CMB temperature anisotropy field T(r), 

dt = ldnr,b{h.pt)T{h) (3-2) 

^ A^^6(r„p,)T,. (3-3) 

j 

where T = Tp is the temperature map in pixel 
basis. Typical mapmakers however ignore the ef- 
fect of extended beam, as it is computationally 
prohibitive, and it "bins" (averages) the observed 
data at each pixel (in the absence of correlated 
noise). So the binned map is made assuming an 
infinitely narrow (Dirac ^-function) beam, in the 
discrete pixel space, which is equivalent of setting 
the beam response to 1 in the pointing pixel and 
outside (that is, a Kronecker S). Then the ob- 



served binned map T = Tp can be written as 



EtAf 



(3-4) 



Here the numerator adds all the detected signal if 
the pointing direction falls in pixel i, and the de- 
nominator counts the number of observations in 
that pixel. However, since the beam is extended, 
the observed map is connected to the true underly- 
ing anisotropy map through the following relation 
[see Eq. (3-3)]: 



Ti :-- 



Et Au dt _ Et Au An Ej K^j.pt) Tj 
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(3-5) 
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(3-6) 



Then, in order to satisfy the definition of effective 
beam, _ 

T =: AQBT, (3-7) 



one must choose 



Bij = 






(3-E 



This is indeed the crude definition of effective 
beam stated in section 2, which directly applies to 
the temperature-only case. In fact, the reason for 
arriving at the expression above in a formal way is 
to aid our understanding for the case of polarized 
detectors, where the definition is not straightfor- 
ward. 

Eq. (3-8) is essentially what we compute for 
temperature only beams, the detailed method of 
how it is done in an efficient way is discussed in 
this section and section 5. 

It is clear from Eq. (3-8) that for each sam- 
ple t, the beam function b must be evaluated at 
all the beam pixels j by rotating the detector in- 
trinsic beam function. We do this by transform- 
ing the pixel center coordinates to the coordinates 
in which the detector intrinsic beam function is 
specified. Note that these operations are the most 
expensive mathematical operations in the effective 
beam computation, as they are used about tril- 
lions times (A^tod x ^beam), SO wc take good care 
in writing them in efficient forms. 



• First, we obtain the coordinates of the point 
in the (locally ffat) orthogonal coordinates 
defined by the polar and azimuthal unit vec- 
tors (^, cj)) with the beam center (r) as the 
origin. 

The radial, polar, and azimuthal unit vec- 
tor in spherical polar coordinates are respec- 
tively given by: 





4> 



smOcos(j)i + sin ^ sin (^j + cos^k 
cos6>cos0i + cos6>sin^j — sin6>k 
— sin^i + cos^j 



(3-9) 



Then, with respect to the beam centre, the 
coordinates (Ax, Ay) of a pixel at r^ can be 
written as: 



Ax 



in -r)e 


(3-10) 


sin6>5 cosO cos((/)5 - 


-0) 


— sin 6 cos 6>5 


(3-11) 


(f5 - r) • 


(3-12) 


sin Ob sin(05 — (j)) . 


(3-13) 



A^ 



We also provide a less accurate but faster 
(by ~ 10% of the whole computation) set of 
formulae to get these coordinates. We have 
tested that these formulae provide enough 
accuracy for practical purposes. 



Ax = 



6>5 C0S(^5 - ^) , 

if < 7r/2 

TT - (tT - Ob) COs((/)5 

otherwise 



(3-14) 



Ay 



smO X 

sin(^5 - ^) , 

ff 105 - 01 > 1° and 

\Ax\ < FWHMbeam 
06-0, 

otherwise . 

(3-15) 

• We then rotate the coordinates by the angle 

'^^ = ip{t) — ippoi + '^eii to transform to the 



beam coordinate system: 



A?/beam / 

COS i/ji) sin i/j}) 
— sin ^|Jl) cos '^5 



Ax 
Ay 



(3-16) 



The beam function is a known input in the 
new Axbeam, A^beam Coordinates. For exam- 
ple, if the intrinsic beam is ehiptical Gaus- 
sian with major and minor axes (Ja,(Ti) re- 
spectively, the value of the rotated beam 
function at 1*5 will be 
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In this paper all the numerical exercises have 
been done assuming elliptical Gaussian beam 
shape with parameters read from a database of 
Planck focal plane parameters and the accurate 
formulae for beam coordinates, Eq. (3-11) and 
Eq. (3-13), have been used. 

3.2. Polarized Effective Beams 

To define effective beams for polarized detectors 
we again formally express the output of a map- 
maker in the noiseless case. This expression is 
then generalized to the extended beam case and 
the polarized effective beam is defined. In case 
the beam is significantly extended, we also need to 
consider the coordinate changes within the beam 
(as the local meridians — the polarization reference 
axes — at different pixels may become significantly 
different). We supply the correction formulae in 
the end of this subsection, though for Planck main 
beams we achieved adequate accuracy without us- 
ing this correction. 

3.2.1. Polarized binned map 

The binning process is independent of the beam 
size. However, if the beam has infinite resolution 
(^-function) the binned map should be identical to 
the true map in the noiseless case. So we derive the 
formula for making binned polarized maps starting 
with such an idealized case, with the view that the 
procedure will be valid for any beam. 

A detector in a CMB experiment records only 
one number for each sample, a linear combination 



of three Stokes parameters /(r), (5(r), ^(r). For 
the idealized detector, one can write 



u;i(Vt)/(r(t)) + 



(3-17) 
wsii^t) U{v{t)) , 



where wt = [wi{^pt),W2itpt),W3{'4>t)] is the weight 
vector at sample t. The above equation can be 
written in a discrete form using the pointing ma- 
trix Afp as 
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For the polarization sensitive instruments used in 
the Planck satellite, it is a good approximation to 
write the weight factors as. 



7 sin(2V^t) 



(3-19) 



where 7 := (1 — e)/(l + e) is the polarization effi- 
ciency, e being the cross-polar leakage. 

Since the detector provides one number for each 
sample, but one needs to estimate three numbers 
(/, Q, U) for each pixel, there must be at least 
three observations per pixel and the set of linear 
equations must be solved. Ignoring the issues of 
degeneracy (not enough independent detector con- 
figurations), below we describe how one can solve 
the equations to make a binned map. 

The goal of making a binned map is that, if the 
beam was indeed vary narrow and there was no 
noise, one should get the true sky map as the re- 
sult. We repeat the procedure presented in jones 
et al. (2007). Multiply Eq. (3-18) through by Wt 
and sum over the samples in a pixel (enumerated 
by the pointing matrix A), to obtain 



^ AtpMVtdt = ^^ AtpAtp'MVtv/J 

t p' \U 

\ / p' 

(3-20) 
Because the pointing matrix has only one unit en- 
try per sample, the sum over p' collapses to: 



^ Atp^tdt = Hp 
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U 



(3-21) 



On the left hand side we have binned weighted 
TOD and on the right we have the true sky times 
the 3x3 polarized "number of observation matrix" 
at pixel p 



H„ 



Y^AtpWtwJ . (3-22) 



To recover the sky, multiply the binned weighted 
TOD by the inverse of number of observation ma- 
trix. So the polarized binned map is given by 

(/ Q uY = H;' ^Atpwtdt 



Yl ^tp^t^^ 



^ AtpV^tdt . 
t 

(3-23) 



As mentioned above, this is the general binning 
formula for any beam size. Also note that, not 
surprisingly, it is similar to the formula for the 
temperature only case [Eq. (3-5)]. 

3.2.2. Polarized effective beams 

We can now calculate the effect of binning on 
a data set observed with a realistic beam of finite 
size and thereby arrive at the formula for polarized 
effective beams. We follow the same path as we 
did for temperature only beams. 

Following Jones et al. (2007), the observed sig- 
nal at sample t is simply the integral of the con- 
traction: 



dt 



Jdnr,b{h.Pt)w^{h.Pt)-{i Q u)l 

(3-24) 
A^^6(n^,p,)W^(n^,p,)-(/ Q /7) J , 



(3-25) 



where W(r5,pt) is the polarization weight vector 
for general extended beam, 

W(r5,pt):=| 7nr6,Pt)cos(2V^,) ), (3-26) 

77^(r5,Pt)sin(2'0t) 

where 7^(1*5, p) is the normalized beam response 
function, defined as 

^(- n\ P\\{h)-P±{h) .^„^. 



expressed in terms of the co- and cross- polar beam 
response functions P||(r6) and P±{rb) respectively. 
7^(1*5, p) can be obtained for any given pointing p 
following the same procedure as that for the beam 
profile 6(1*5, p) described explicitly for temperature 
only beams. 

On substitution into the mapmaking equation 
[Eq. (3-23)], the observed (convolved) sky map be- 



H-^ ^A,,w, X (3-28) 



Al]^6(r„p,)W^(n„p,)- Q 

Thus, in order to satisfy the definition of polarized 
effective beam, 

(7 Q U)J =: Al^^Bi,-(/ Q U)J, 

3 

(3-29) 
one must choose, 




Br 



^ AtpWtwf 



(3-30) 



^ Ati b{rj,pt) WtW^{hj,pt) . 
t 

The effective beam for each pointing pixel i 
now has two parts. There is a 3 x 3 symmet- 
ric matrix H^, which is a generalized version of 
number of observation at that pixel. Then, for 
each neighboring pixel j there is a 3 x 3 ma- 
trix ^^A^i6(r^-,pt)wtW^(n^-,pt), which encap- 
sulates the beam profile and need not be a sym- 
metric matrix in general. Computation of this 
quantity can be done in a similar way as in the 
temperature only case, first computing the coor- 
dinates of neighboring pixels with respect to the 
beam center and then evaluating the polarized 
beam profiles 6(rj,pt)W^(nj, p^) using the mea- 
sured/modeled intrinsic detector response func- 
tions. 

The cross polar beam is normally much smaller 
than the co-polar beam and, in practice, the co- 
polar beam might be the only response which 
could be accurately measured. These facts lead 
to a reasonable assumption that Virh) ^ I7 so 
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that, W(r5,pt) = Wt, which, in turn, simphfies 
the computation of effective beams. 
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(3-31) 

as the beam profile for all the polarization becomes 
the same, except for weight factors (not depen- 
dent on f-j), and makes all the component matri- 
ces symmetric. Note that, if the effective beam 
matrices are not symmetric'^, computational stor- 
age requirement increases by 50% (9 instead of 6 
components for each matrix). Computation time 
is not much affected by this modification. 

FEBeCoP does include the ability to compute 
effective beam matrices of both kinds, symmetric 
or otherwise. However, the numerical results and 
comparisons presented in this paper are for the 
possibly more practically achievable scenario of 
V{Yh) = 1- So, unless otherwise specified, the ma- 
trices representing polarized effective beams com- 
puted for accuracy and performance analysis of 
this method would be symmetric. 

3.2.3. Transformation from beam frame to sky 
frame 

If the beam is significantly extended, the pixels 
which are far (that is on very different longitudes) 
would suffer from an extra polarization rotation 
effect, this is because the far away longitudes are 
no longer close to parallel (except when they are 
very close to the equator) , so the frame of polar- 
ization also rotates. The mathematics to account 
for that extra rotation is presented below. 

If the detector orientation at sample t is repre- 
sented by Euler rotation matrix 'D{(j)t^6t,i^t)^ then 
the polarization field for the beam in pixel j is 
evaluated at a direction D~^(0t7 ^t^ i^t) • n^ in the 
beam frame. The orthonormal polarization bases 
differ between frames by a rotation. Stoke / is a 
scalar, unchanged by the rotation. To transform 
(Q, U)^ we need to know the angle between the po- 
larization bases in the sky and beam frame (Fig- 
ure 3). Imagining the spherical triangle between 



^ "Polarized effective beam matrices are symmetric" signifies 
that the 3x3 matrix for representing polarized effective 
beams is a symmetric matrix, it should not be confused 
with an effective beam constructed from symmetric detec- 
tor beams. 
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Fig. 2. — Focal plane layout for Planck (Dupac 
.^ Tauber 2005; Tauber et al. 2010), showing the 
position, approximate shape, and orientation of 
the beams of every detector. Crosses show the 
orientation of the polarization axes of the polar- 
ization sensitive detector pairs. The arrow on the 
left shows the scan direction, which is fixed for the 
duration of the mission. 
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Fig. 3. — Orientation of the detector for sample 
i after rotation by Euler angles {0^,^^,'^^}, and 
the orthonormal basis (0, 6) at pixel p in both the 
sky frame (where 6 points along the great circle to 
the pole ) and the beam frame (where 9 points to 
the detector). The angle between the polarization 
bases is labeled Cip- 
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the North pole, the position of sample i, and pixel 
p, we may deduce the angle between the orthonor- 
mal bases at the pixel location, denoted Cip, with 
two applications of the law of haversines, 

haversin(A$^p) = haversin(6>^ — Op) + 

sm{6i) sin ( ^p) haver sin (^^ — (j)p) (3-32) 

haversin(^i) = haversin(^p — A^ip) + 
sin(^p) sin(A<l>^p)haversin(C^p), (3-33) 

where haversin(^) = sin^(^/2). The first equation 
solves for A^ip, representing the great circle dis- 
tance from the pixel to the sample position, which 
is required in the second equations to solve for dp. 
This extra rotation would lead to a transforma- 
tion of the observed sky in the following way. 



u' 



cos2Ctp sin2Ctp f Q 
-sm2Ctp cos2Ctp \ V U 

(3-34) 
This would change the beam weight factor in the 
expression for effective beams to W(np,pt) -^ 
W'(np,p^) = W('0t + Ctp^Pt) and the expression 
would become 
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(3-35) 



^ Ati b{rj,pt) w^W'(n_^-, pt 



It is easy to see that the effective beam matrix at 
each pixel is no longer symmetric in this case, even 
if the cross-polar beam is ignored. 

4. Theoretical Accuracy 

Effective beams defined over a grid could have 
quadrature accuracy issues if the pixelization was 
not fine enough. Also, the radial cut off we im- 
pose on the beam leads to some information loss. 
In this section we study the truncation error in in- 
tegrating a 2-D Gaussian over a HEALPix grid 
to determine the minimum value of the cut-off 
radius and also show that quadrature error is 
sub-dominant here. Then, by comparing effective 
beam convolved CMB maps using Total Convolu- 
tion for a toy scan model, we show that the default 
map pixel resolution (discussed in section 2) pro- 
vides sufficient accuracy. 



4.1. Effective Area Loss of a Truncated 
Gaussian 

By definition, the surface integral of a 2-D 
Gaussian over an infinite plane is unity. How- 
ever, if we cut off a zero mean Gaussian with 
variance a^ at rmax = ^fwhm x FWHM = 
2nFWHM a/2 In 2 a, the "area loss" turns out to be 
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e := 1- 
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The order of magnitude error is thus 
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(4-2) 
That is, if we choose the cutoff at npwHM = 2, 
the area loss will be ~ 10~^, which should pro- 
vide enough accuracy for convolution. It is worth 
noting that in the context of CMB maps the error 
is actually the difference between convolved val- 
ues with and without cutoff, which is dependent 
on the underlying CMB map. Since the maps have 
zero expectation, the expected error would also be 
zero, but the fractional RMS error would be the 
same as the effective area loss. 

In practice convolution is done in discrete space. 
To confirm that the above analytical estimate still 
holds and the usual choice of pixel resolution does 
not add any significant quadrature error, we do the 
exercise numerically. We compute the area loss for 
a Gaussian with FWHM - 32' (Planck LFI 30GHz 
beam size) using two sets of beam cut-off radius, 
^FWHM = 2 and npwHM = 4, and pixelization, 
A^pix = 12 X 512^ and TVpix = 12 x 2048^ 

Figure 4 shows plots of area loss for the four 
cases as functions of latitude (the evaluation was 
done along the meridian without any loss of gen- 
erality). Clearly, area loss depends on the cut- 
off radius and can not be improved by increasing 
A^pix to higher than the usual value for the given 
beam size. There is some quadrature error near 
the poles, but we can ignore that as the solid an- 
gle which is affected by this is very small (only 
few pixels). There is a sharp features near each 
HEALPix polar circle for higher cutoff radius, only 
if the target accuracy is better than 10~^. Our de- 
fault choice would be npwHM = 2, which provides 
adequate accuracy without bringing in any signif- 
icant quadrature issues. 
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4.2. Comparison with Total Convolution 
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Fig. 4. — A symmetric Gaussian beam of 
FWHM=32.2352' was integrated in the pixel space 
for two different HEALPix pixel resolution as well 
as two different radial cut-offs (expressed in the 
units of FWHM of the Gaussian). In the ideal 
case the result of the integration — the "area" of 
the Gaussian — should be unity. The deviation of 
the numerical results from unity, the "area loss", 
is shown in the plot. As expected, the area loss is 
mainly dependent on the cut-off radius and very 
weakly dependent on the quadrature error due to 
pixelization. The sharp features in the plots are re- 
lated to the pixelization scheme (HEALPix). The 
point to note here is that, for a cut-off radius of 
^FWHM = 2 and ngide = 512 the errors are at 
the level of single precision floating point accu- 
racy, which we consider sufficient for the target 
accuracy of map simulation. 



Next we compare effective beam convolved 
maps for toy models of beam orientations with To- 
tal Convolution Data Cubes for different choices 
of pixel resolution. TCDC are theoretically exact 
for band limited beam functions on an equi-spaced 
grid of beam location and orientation. We used 
those values of orientation in the effective beam 
computation and convolved with the same map (in 
pixel basis) from which the TCDC were derived. 
TCDC are computed at equi-spaced intervals of 
RA and dec, so the beam centers do not coincide 
with the HEALPix pixel centers. This is similar to 
the real scans, which need not pass right through 
the pixel centers. One small approximation is 
unavoidable here, a Gaussian beam is not purely 
band limited, specially for asymmetric beam there 
is a cut-off on the azimuthal index ^m\ This error 
is not severe unless the accuracy requirement is 
very stringent. In fact, this error adds to true dif- 
ferences making the test even more strict for the 
pixel space convolution. 

Figure 5 shows that for a ~ 32' beam, the nat- 
ural choice of pixel resolution ngide = 512 pro- 
vides very good accuracy and further increase in 
resolution does not improve it by any significant 
amount. At low ngide the beam is very poorly 
sampled, which leads to orders of magnitude bet- 
ter accuracy for each step of ngide improvement. 
However, beyond ngide = 512 the amount of ex- 
tra information added is low, so the accuracy im- 
provement is marginal and is also limited by single 
precision floating point accuracy. So, we may con- 
clude that the usual choice of pixelization is the 
optimum choice. 

5. Computation Algorithm and Imple- 
mentation 

In the previous sections we have presented the 
theoretical statement of the problem and error 
analysis. The practical implementation of the 
method is a bigger challenge due to the organi- 
zational complexities of efficient distributed com- 
putation. This section outlines the possible routes 
of computation, their relative merits and demerits 
and some final outcomes. 
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5.1. Quantities to Compute 



Difference between TotalConvolver and pixel space convolved maps 
for Planck LFI-28a beam, cut-off radius 2 x FWHM 
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Fig. 5. — Difference between Total Convolution 
maps (data cube slices) and pixel space convolved 
maps for different HEALPix ngide is shown in 
this figure for symmetric (dashed) and asymmet- 
ric (solid) beams of FWHM - 32' and elliptic- 
ity ~ 1.4. Total Convolution is exact on an equi- 
spaced (6>, (j)) grid for band limited functions. The 
comparisons were performed on this grid to access 
how far pixel space convolution is from the exact 
theoretical convolution and to determine the op- 
timal ngide foi" effective beam computation. For 
each HEALPix pixel, the TCDC value at the clos- 
est grid point was used for comparison, no inter- 
polation was performed. 



The definitions provided in section 2 suggest 
that in order to evaluate and store the effective 
beams, the following quantities have to be com- 
puted for each pixel: 

• The list of pixels p that are within the beam 
cut-off radius from the pixel center. Ideally 
one should find this list with respect to the 
sample beam centers within the pixel, but 
the difference would be negligibly small and 
cost large amount of computation, as it will 
have to be done at each sample. Instead, the 
list is computed in the very beginning with 
respect to the pixel centers. The number of 
element in the list, A/'ngt, is slightly different 
for each pointing pixel. 

• The number of observation A^obs at a given 
pixel or, for polarized beams, the weight ma- 
trix H (whose first component is A'obs)- Note 
that effective beams are computed in a way 
that they are always normalized (divided 
by A/obs) as they are developed. That is, 
only the fractional changes from the existing 
beams are added. This prevents large num- 
ber arithmetic (leading to numerical errors) 
near the ecliptic poles where, for Planck, 
Aobs is higher by 3-4 orders of magnitude. 

• The effective beam Bij at the neighboring 
pixels for temperature only beam and the 
effective beam matrix Bij for each neighbor- 
ing pixel for polarized beams. 

In addition, there can be intermediate quan- 
tities that need to be pre-computed for conve- 
nient execution of the main steps, this intermedi- 
ate quantities can also consume significant amount 
of computation. Details of these implementation 
dependent steps are provided in the following sec- 
tions. 



5.2. Algorithm 

From the description provided in section 2 it 
follows that there are two obvious ways to com- 
pute the effective beams, using a time-domain or 
pixel-domain parallelization scheme. Their differ- 
ence is purely computational, the challenge being 
to manage the large effective beam data object and 
divide computation efficiently. We explored both 
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approaches while developing our effective beam 
pipeline, settling on the pixel domain as our pre- 
ferred method of parallelization. 

5.2.1. Time domain parallelization 

Time domain parallelization means different 
portions of the scan are assigned to different pro- 
cesses. As the scan progresses, more and more 
pixels are assigned to each process and memory 
is allocated progressively. To avoid memory over- 
flow, after certain predetermined interval, beams 
are written to disk — if for some pixels beams al- 
ready existed on disk, they are read, new beams 
are added and rewritten to disk. 

• Advantages: Simplistic implementation 
(though the memory management process 
is complicated) 

• Disadvantages: I/O intensive in different 
ways. Disk I/O is many times the total size 
of beams, as partial beams are read and writ- 
ten many times. Moreover, since one pixel 
can appear in different processes, the total 
beam (data) size is about three times bigger 
than the true size for one detector, leading 
to a massive loss of memory, disk space and 
overall I/O overhead; handling multiple de- 
tectors at a time is practically impossible, 
as rings would become "thicker" due to the 
angular distance between detectors, requir- 
ing much more memory for each ring; it is 
not easy to locate a pixel in the beam files 
which makes it very difficult to load beams 
for a certain set of pixels in memory - hence 
the beam files have to be read from disk 
for every convolution; there is no easy way 
to predetermine memory usage by a certain 
process, so memory efficient parallelization 
is practically impossible. Finally, each pro- 
cess requires the satellite pointing informa- 
tion simultaneously, which in turn requires 
reading the same satellite pointing file — this 
is a big bottleneck and increasing the num- 
ber of processes does not necessarily speed 
up computation. 

Due to the long list of disadvantages above, this 
method was not very useful for efficient compu- 
tation of effective beams, especially for polariza- 
tion. Though the simplicity of the implementation 



helped us making initial assessment of accuracy, 
computational requirement and feasibility of the 
method. 

5.2.2. Pixel domain parallelization 

Pixel domain parallelization means different 
pixels are assigned to different processes, though 
the division of pixels need not be equal as we ex- 
plain below. 

Since the satellite scan path is efficiently ob- 
tained in the time domain with contiguous stretch 
of samples, we need to go through an intermedi- 
ate step of binning the detector angles in pixels — 
for each pixel we build a list of three scan an- 
gles (6>, 0, ip) for all samples in that pixel, in other 
words, the whole scan pattern is transformed to a 
set of "Pixel Ordered Detector Angles" (PODA). 
Though this step requires significant amount of 
computation (few x 10% of total computation), 
once the satellite pointing information is available, 
this binning is fixed, so this step need not be re- 
peated. This means the effective beam for dif- 
ferent beam shapes can be recomputed without 
repeating this step. Also, polarized detectors in 
Planck are in pairs whose pointing angles are the 
same and polarization angle differ by a constant 
threshold, so one need not compute PODA for half 
of the polarized detectors. 

The next step is the load division among com- 
puting processes. Though assigning equal number 
of pixels to each process would lead to the most 
efficient division of memory, we assign equal num- 
bers of time samples to each process for the fol- 
lowing reason: normally the experimental scans 
are designed in such a way that nearly the whole 
sky is uniformly sampled, which ensures that even 
distribution of samples would distribute most of 
the pixels evenly. However, some pixels (pixels 
near the ecliptic poles for Planck) can have orders 
of magnitude more observations than the average. 
These pixels would slow down the corresponding 
processes by orders of magnitude and, as a result, 
the whole computation. Even division of pixels 
along longitudinal strips would not help, as the 
distribution of these pixels is highly asymmetric. 
But in the equal number of samples scheme, these 
few pixels will be assigned to processes given many 
fewer pixels than average. This does not signifi- 
cantly increase the number of pixels in the other 
processes, since the number of these highly ob- 
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served pixels is low. 

• Advantages: Small I/O overhead, the 
whole beam is kept loaded in memory, and 
optionally written to disk only in the end; 
each process reads and buffers POD A sep- 
arately and no bottleneck situation occurs; 
the process is highly scalable, can use thou- 
sands of processors; memory management is 
transparent and possible to keep the whole 
beam for multiple detectors loaded in mem- 
ory for the convolutions, without ever writ- 
ing to the disk. 

• Disadvantages: Complicated implementa- 
tion; requires significant preprocessing and 
disk space to store PODA. 

This method is highly efficient and has been 
implemented for Planck mission. From this point, 
our effective beam computation refers to pixel 
space parallelization, unless mentioned otherwise. 

5.3. Implementation for Planck 

We have fully implemented the effective beam 
method for the Planck mission. The mission 
specific inputs required for the effective beam 
computation are the detector characteristics and 
the pointing information. Planck's detector 
characteristics — beam parameters, polarization 
angles, cross polar leakage etc. — are read from 
a focal plane database. The detector pointing 
information is extracted from compressed satel- 
lite pointing using the Generalized Compressed 
Pointing (GCP) library through the MADCAPS 
interface ^. 

The beams are stored in binary formats which 
are accessible to both C and Fortran standard I/O 
routines. Each process writes its own binary file, 
effective beams and other data for each pixel are 
stored as one "record" in one of those files. To 
preserve compatibility between different comput- 
ing languages, the record size is kept the same for 
each pixel, which is determined by maximum A^nst- 
The address of a pixel is stored in a combined in- 
dex file, which is actually a sky map — each map 
value encodes the file number and location for that 
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pixel as a single precision integer. A software in- 
terface to conveniently load partial or full effective 
beams has also been developed. 

The cost of effective beam computation scales 
as A^tod X ^beam, whcrc A^tod IS the total number 
of samples (sum of number of observations at at 
ah pixels). For 15 months of Planck HFI 143GHz 
channel (simulated) data, sampled at 172Hz (to- 
tal 6.8 X 10^ samples), with TVpix = 12 x 2048^ 
pixels, the computation cost is ~ 130 CPU hours 
per detector. For Planck LFI 30GHz channel, the 
cost time is ~ 30 CPU hours. These costs could 
be reduced by ~ 10% by using the little less accu- 
rate formulae for beam computation provided in 
section 3.1. 

The total size of the effective beam can be es- 
timated as follows: for each pixel with Nnst beam 
pixels, total volume of data is essentially the num- 
ber of neighboring pixels plus number of beam ma- 
trix elements. So for polarized beams, if we use 
single precision numbers, the size of each record 
is approximately 4(1 + 6)A'iist = 287Viist Bytes, so 
the total size of the polarized effective beam will 
be 28A^iistA^pix Bytes. For Planck HFI 143GHz 
with nside = 2048 and N^st ^ 225 the total vol- 
ume of effective beam turns out to be ~ 290GB. 
This is approximately the disk space required to 
store the beams, the total amount of (distributed) 
memory requirement is greater than this, but of 
the same order. Note that, if we chose to use dif- 
ferent beams for unpolarized and polarized com- 
ponents, the polarized effective beam matrices will 
not be symmetric, and 50% more storage would be 
required. The other factor which may affect the 
size is the cut-off radius. Increasing this may im- 
prove the accuracy slightly for polarized detectors, 
but the total beam object size scales as the square 
of the radius, so the increase in this data volume 
can be significant. 

We illustrate the features of effective beams and 
PSFs derived from them by showing a few figures 
for Planck 30GHz (simulated data). The distri- 
bution in the number of observations due to the 
scan pattern is shown in Figure 6. In Figure 7, 
we show logarithmic plots of effective beams (left) 
and PSFs (right) at different locations on the sky, 
where the locations have been marked in Figure 6 
(in Galactic coordinates) respectively by "cross" 
(ecliptic North pole) , "square" (a cusp near ecliptic 
North pole), "triangle" (a mid-declination point) 
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and "asterix" (ecliptic equatorial point). We over- 
lay the isocontours at 50%, 10%, 1% and 0.1% 
levels for the numerical maps and also elliptical 
Gaussian fits to the beams and PSFs. The fitted 
values, FWHM, ellipticity, orientation angle with 
respect to the local meridian (vertical axis) and 
the normalized RMS fit error have been mentioned 
in the subtitle. The center of the central pixel (ef- 
fective beam/PSF center) and the center of the 
fitted Gaussian have been marked with "+" and 
"x" respectively, when they coincide the markers 
form an "eight arm asterix". These figures tell 
us that the effective beam / PSF has significant 
variations across the sky, implying that simplify- 
ing assumptions for asymmetric beam studies may 
not be very accurate. Specially, near the ecliptic 
poles and the cusps of the scans, the asymmetry 
can be large. 

6. Application of Effective Beams - I: Con- 
volution 

The main purpose of the effective beams is very 
fast convolution and the FEBeCoP implementation 
handles it efficiently. Below first we provide some 
approximate computation cost estimates for the 
effective beam convolutions; next we test the ac- 
curacy by comparing with the results derived using 
existing Planck simulation software "level-S" and 
analyzed by HEALPix utility anafast (Gorski 
et al. 2005). 

6.1. Efficiency 

Convolution with effective beams not only 
scales as A^pix x A'beam, the proportionality con- 
stant is also very small — of order 1 for temperature 
only beam and of order 10 for polarized beams. So 
the convolution is very fast, but it is still impor- 
tant to ensure that the implementation is efficient 
and free of I/O bottlenecks. Among other opti- 
mizations, we pre-invert the polarization "hits" 
matrices and ensure that the beams are read in 
contiguous chunks in order to minimize file ac- 
cess. Our current implementation is now solely 
limited by map I/O — this is indeed satisfactory 
that the core computation is faster than the time 
taken to read a map from disk, but, on the other 
hand, unfortunate, that we can not make it faster! 
However, studies that involve online simulation 
and analysis of maps, e.g., transfer function es- 



timation, where maps need not be read from or 
written to disk, run almost at the native speed, 
though they involve some amount of network I/O. 
We provide some typical computation costs in ta- 
ble 1 for different beam sizes, assuming that the 
usual number of HEALPix pixels were used to 
represent a map. 

6.2. Accuracy: Comparison w^ith Existing 
Planck Simulations 

Next, we assess the accuracy of the convolution 
by comparing with the simulated maps generated 
using the standard Planck Level-S TOD simula- 
tion pipeline in combination with the Springtide 
or MADAM map maker (Ashdow±± ^ ). 

Starting from the focal plane database, the satel- 
lite scan strategy and (the spherical harmonic 
transforms of) an input map, we generate To- 
tal Convolution Data Cube with r?imax = 2 for 
symmetric beams and mmax = 14 for asymmet- 
ric beams. The TCDC are then interpolated by 
multimod (Reinecke et al. 200(), to generate TOD 
at each scan sample. These (noiseless) TODs are 
then binned by a mapmaker to generate maps. 
The effective beams are computed using the same 
pointing and focal plane database and convolved 
with (the pixelized version of) the input map. The 
comparison between the maps and the power spec- 
tra for two Planck channels, 30GHz (all detec- 
tors: 27a, 27b, 28a, 28b) and 143GHz (four de- 
tectors: la, lb, 4a, 4b), are shown below. Among 
all the Planck beams, the 30GHz ones have the 
highest asymmetry (ellipticity ~ 1.35 - 1.40) and 
the 143GHz ones have the lowest asymmetry (el- 
lipticity - 1.02 - 1.11). 

We first show simulated full sky temperature 
and polarization maps for Planck 30GHz channel 
using Level-S and FEBeCoP and their differences in 
left and right panels of Figure 8. The maps from 
the two methods are clearly similar and the differ- 
ence is very small — the fractional RMS difference 
for temperature maps is few x 0.01% and for po- 
larization amplitude the difference is few x 0.1% 
(note that the contrast in the difference maps 
have been amplified by using a histogram equal- 
ized color map). We have also plotted zoomed 
versions of the above maps for North polar and 
equatorial regions in Figure 9. The agreement 
between FEBeCoP and levelS is good in all the 
places, expect for only few pixels (among millions) 
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Beam 
FWHM 

(arcmin) 


HEALPix 

^side 


Temperature only 


Polarized 


Effective 
beam size 

(GB) 


Computation time 


Effective 
beam size 

(GB) 


Computation time 


beam/detector 
(CPU hr) 


convolution 
(CPU min) 


beam/detector 
(CPU hr) 


convolution 
(CPU min) 


32 


512 


7 


27 


1.33 


24 


31 


4.4 


13 


1024 


18 


42 


5 


64 


48 


17 


7 


2048 


90 


128 


21 


293 


132 


64 



30GHz § observation 



Table 1: Approximate computational resource estimates for FEBeCoP to compute effective beams for ~ 15 
months' of observation (with cut-off radius of 2 x FWHM) and to convolve using the (precomputed) effective 
beams on the NERSC computing facility. 



right at the poles, which is not significant. (We 
omit these plots for 143 GHz, where they are un- 
informative.) The power spectra of the maps from 
two methods for both 30GHz and 143GHz^ chan- 
nels are shown in [Figure 10 and 12], the spectra 
are almost indistinguishable (the lines nearly co- 
incide). Their fractional differences reveal more 
details [Figure 11 and 13]. What is important to 
notice in these comparisons is that the difference 
between the effective beam convolved maps and 
the Level- S plus mapmaker maps, and their power 
spectra, is much smaller than the corresponding 
difference between a Level- S and symmetric Gaus- 
sian smoothed map, obtained using the HEALPix 
utility synfast (Gorski et al. 2005). Especially 
for polarization the symmetric beam assumption 
systematic is large (~ few 10%), while the effec- 
tive beam error is < 0.1%. Even for symmetric 
detector intrinsic beams, small asymmetry is in- 
troduced to the effective beams by scan and that 
too is faithfully reproduced by FEBeCoP. This can 
be seen in the top panel of Figure 10. Here the 
spectrum of theoretical symmetric Gaussian con- 
volved map computed using synfast differs from 
the spectrum of levelS map with symmetric beam 
and simulated Planck scan, but the spectrum of 
FEBeCoP map for the same inputs lies exactly 
on top of the LevelS spectrum. These exercises 
demonstrate that effective beam convolution pro- 
vides very good accuracy in mimicking the sys- 
tematic effects in the detectors. 




Log 



Fig. 6. — Log plot of typical number of observa- 
tions at each pixel of simulated Planck sky maps 
for ~ 15 months' scanning. The one shown in 
the figure is the total number of observations for 
four 30GHz detectors in Galactic coordinates. The 
marked points correspond to the locations of beam 
and PSF plots shown below. 



^We observed that for the 143GHz channel, increasing the 
cut-off radius from 2 times FWHM to 2.25 times FWHM, 
increases the polarization accuracy by a small amount. 
However, the storage also increases by 25%, so there is 
some small trade-off involved in the choice of parameters. 
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30GHz eff beam: pix #0 30GHz PSF: pix ^0 
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(46,0, ae.e) Ecliptic 



(46,0, ae.e) Ecliptic 



30GHz eff beam: pix #52836 3QGHz PSF: pix #52836 
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30GHz eff beam: pix #523264 
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30GHz PSF: pix #523264 

Fit:ep„5i,=3a.5a',e=1.317,f,„=-4.S*,A„ijgB,= l-'^l% 
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30GHz eff beam; pix #1572664 30GHz PSF; pix #1572864 

Fit:ej,,^=32.56',e=l,3S5,f^=0,e°,VgE=0-36% Fit:er,^=32.52',e=l,3S6,f^=0,7°,VgE=0-6! 




levels + madam: T 



levels + madam: \Q,U\ 




Fig. 8. — Comparison between simulated temper- 
ature (left panel) and polarization (right panel) 
anisotropy maps for Planck 30GHz channels ob- 
tained using Level-S (top) and FEBeCop (middle). 
The maps are visually quite similar, the difference 
map (bottom) are also shown using a histogram 
equalized color map (to amplify the contrast in the 
image). The fractional RMS difference for temper- 
ature is 0.04% and for polarization is 0.9%. 



(iao.1., 0.0) Ecliptic 



Fig. 7. — Log plots of effective beams (left) 
and PSFs (right) at four different locations of 
sky - north ecliptic pole, a cusp near North 
ecliptic pole, a mid-declination point and eclip- 
tic equator, marked in Figure 6 respectively by 
"cross" , "square" , "triangle" and "asterix". The 
isocontours of the map and elliptical Gaussian fit 
to the map are overlayed at 50%, 10%, 1% and 
0.1% levels respectively. The center of the effec- 
tive beam/PSF and that of the Gaussian fit have 
been marked respectively by "+" and "x". 
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Fractional difference between Angular Power Spectra for LFI 30GHz 
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Fig. 10.— Comparison of TT (top), EE (middle) 
and TE (bottom) power spectra for LFI 30GHz: 
the main purpose of this plot is to show that Level- 
S and FEBeCoP spectra match very well for both 
both symmetric and asymmetric beams (the lines 
are on top of each other) . This plot shows that if 
the detector intrinsic beam is perfectly symmetric, 
the spectra are close to those of synfast generated 
symmetric Gaussian smoothed maps. 
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Fig. 11. — Difference between power spectra of 
the maps and power spectra of difference maps for 
LFI 30GHz. While the difference between Level-S 
maps and ideal symmetric beam convolved map 
using synfast is large (dashed), few 10%, the dif- 
ference level between FEBeCoP and Total Convolu- 
tion + Level-S + MADAM is less than -- 0.5% at 
almost all multipoles (solid). This shows that the 
effect of asymmetric beam for polarization is large 
and FEBeCoP can imitate the effect accurately. 
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Comparison of Angular Power Spectra for LFI 143GHz 
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Fig. 12.— Comparison of TT (top), EE (middle) 
and TE (bottom) power spectra for HFI 143GHz: 
the main purpose of this plot is to show that Level- 
S and FEBeCoP spectra for both symmetric and 
asymmetric match very well (the lines are on top 
of each other). This plot also shows that if the 
detector intrinsic beam is perfectly symmetric, the 
spectra are close to those of synfast generated 
symmetric Gaussian smoothed maps. 



Fractional difference between Angular Power Spectra for LFI 143 GHz 



levels - synfast 
Cj of (levels - synfast) map 
- FEBeCoP - levels 
Cj of (FEBeCoP - levelS) map 




200 400 600 800 1000 1200 1400 1600 1800 2000 

1 



Fractional difference between Angular Power Spectra for LFI 143 GHz 



0.25 

0.2 
0.15 " 

0.1 



levels - synfast 
C, of (levels - synfast) map 
- FEBeCoP - levels 
Cj of (FEBeCoP - levelS) map 



A A 



^^^t/n^ 



/\ 



'' ,..-i 



trt 



'/ 



\£ 



200 400 600 800 1000 1200 1400 1600 1800 2000 

1 

Fig. 13. — Difference between power spectra of 
the maps and power spectra of difference maps for 
HFI 143GHz. While the difference between Level- 
S maps and ideal symmetric beam convolved map 
using synfast is large (dashed), as much as 20%, 
the difference level between FEBeCoP and Total 
Convolution + Level-S + MADAM is less than 
^ 0.02% at almost all multipoles (solid). This 
shows that the effect of asymmetric beam for po- 
larization is large and FEBeCoP can imitate the 
effect accurately. 
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7. Application of Effective Beams - II: To 
Study Collective Effects of Beams 

Effective beams are not merely useful for fast 
convolution, they can provide a more realistic way 
of understanding and interpreting systematic ef- 
fects. To demonstrate the potential applications, 
we provide two examples in this section. 

7.1. Statistics of the Effective Beams 

Complex combination of scanning, beam and 
pixel shapes lead to different shapes of effective 
beams and PSFs, few examples were shown in 
Figure 7. Software tools developed for FEBeCoP 
readily allow one to obtain useful statistics of the 
beams and PSFs. In Figure 14 we show his- 
tograms characterizing the variation of elliptical 
Gaussian fit parameters to the 30GHz simulated 
effective beams, sampled at HEALPix pixel cen- 
ters of nside = 16 (this method ensures uniform 
sampling, though it has a very high probability 
of missing the "cusp" regions in the scan pattern, 
where the beams/PSFs can be drastically differ- 
ent). Typically the effective beams are few per- 
cent bigger than the intrinsic detector beams as 
their centers are spread over a pixel, adding extra 
smoothing. 

7.2. Understanding Observed Power Spec- 
tra 

In the previous section we presented compar- 
isons between symmetric and asymmetric beam 
convolved maps. The observed power spectra of 
these maps were different by few percent. There 
is, however, no trivial way of interpreting this 
difference — the reason is that the effective beams 
have varying shapes across the sky. Previous at- 
tempts to explain this systematic deviation us- 
ing mathematical models assume constant beam 
shape (with possibly rotation) across the sky, 
though the real scan strategy does not have a sim- 
ple analytical model. This problem can be par- 
tially or completely tackled by studying the ef- 
fective beams. Here we use a leading order ap- 
proximation to show that the ratio of tempera- 
ture power spectra observed using asymmetric and 
symmetric beams can be modeled using the ro- 
tated spherical harmonic transforms of the effec- 
tive beams. 

For a symmetric beam of spatially non- varying 
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Fig. 14. — Variation of elliptical Gaussian fit pa- 
rameters to 30GHz simulated effective beams is 
shown in this figure. FWHM is in arcmin and the 
orientation angle is in degrees. The average ori- 
entation angle is close to zero in this plot due to 
the overall detector geometry specific to this case. 
Tools developed for FEBeCoP allows fast and easy 
applications to study the statistical properties of 
effective beams. 
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size, it is trivial to show that the observed power 
spectrum Ci is related to the true power spectra 
through the relation 



c. 



Bi Ci 



(7-1) 



where Bi are the Legendre transforms of the beam. 
In general, though, the effective beam shapes are 
not constant and they vary across the sky, as 
shown in Figure 7. So the transfer function is no 
longer as simple as Bf. The easiest way to esti- 
mate the transfer function is by performing Monte 
Carlo simulations and estimate the ratio Ci/Ci for 
each multipole, which is, in fact, one of the moti- 
vations for the effective beam approach. However, 
to extract more insights directly form the effective 
beams, we crudely estimate average Bf of effec- 
tive beams by taking some sample effective beams 
from different latitude, computing the rotated har- 
monic transform and incorporating the leading or- 
der correction following the algebra in Mitra et al. 
(zui ). 

Even though these approximations are crude, 
the results look quite satisfactory. In Figure 15, 
we plotted the ratio of power spectra of simulated 
maps for Planck 143-5 channel with asymmetric 
and symmetric beams. For the same two cases, the 
ratio of Bf estimated from effective beams is over- 
laid on the plot. Clearly the match is quite good 
for demonstration purposes. There are accuracy 
issues at high multipoles due to many reasons — 
for example, the formula we used is valid for 
beams with reflection symmetry, though the effec- 
tive beams do not have that symmetry — the effec- 
tive beams are constructed from elliptic or circular 
beams which have that symmetry, but since the 
beam centers and orientations are not distributed 
uniformly in a pixel, the reflection symmetry is 
slightly broken, which may introduce extra first 
order correction terms. 

The ratio of B^ of effective beams at certain 
latitudes (not averaged) are also overlaid to show 
that the characteristics of the effective beams vary 
significantly at different sky locations, which may 
break simplified assumptions often used in simu- 
lations. The variation of B^ with latitude need 
not represent any trend, as the effective beams 
can vary significantly from pixels to pixels and 
also planck scanning is not completely azimuthally 
symmetric. 



8. Conclusions 

Realistic simulations of high precision experi- 
ments play an extremely important role in study- 
ing systematic effects of the instrument. There is 
no exception for scanning CMB experiments, like 
Planck. In the last several years there has been 
significant success in developing an end-to-end 
(E2E) simulation pipeline to simulate very realis- 
tic Planck time-order-data (TOD) going through 
the full scan strategy. However, such simulations 
involve a very detailed computational procedure, 
and are usually computationally prohibitive in ap- 
plications to adequate Monte Carlo studies (that is 
repetitive many thousand runs, or more). Many 
CMB analyses do not require full TOD for the 
study of systematic effects, as simulated realistic 
observed maps provide all the ingredients required 
for these studies. In fact, these studies may require 
many more simulations than it is feasible to ob- 
tain by performing the E2E simulations, and then 
going through mapmaking. Beam asymmetry is 
one of the most important potential instrumen- 
tal systematic effect in Planck. Pseudo-C/, and 
other power spectrum estimators could make use 
of large number of realistic sky map simulations 
to estimate the correction needed to estimate the 
true power spectra. This was the motivation for 
developing a pixel based convolution method using 
effective beams, which we have implemented using 
the computer code FEBeCoP. The method precom- 
putes the effective beams and makes it possible to 
perform very fast convolutions to simulate realis- 
tic sky maps. FEBeCoP also provides point spread 
functions (PSFs) at every pixel on the sky, which is 
otherwise non-trivial to obtain. We have success- 
fully constrained the computation cost to a mod- 
erate value as well as managed to properly utilize 
the available distributed computing power. The 
computation requirement for FEBeCoP is modest 
as compared to the modern scientific computing 
facilities available for Planck analyses. Due to the 
simplicity of FEBeCoP interface and its efficiency, 
we believe that it would be very easy and useful 
to incorporate effective beam convolution in many 
systematic effect studies. 
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Comparison of the ratios of pseudo-Cj-s and the effective beam Bj factors 
computed for the symmetric and asymmetric intrinsic detector beams 

Ratio of pseudo Cj from Level-S _ 
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Fig. 15. — Illustration of approximate model- 
ing of the transfer function using harmonic trans- 
forms of the effective beams. Dotted: the ratio 
of observed power spectra of simulated maps for 
Planck 143GHz channel #5 for asymmetric and 
symmetric detector intrinsic beams. Solid: Ratio 
of average effective B^ obtained from the effective 
beams for the same simulations, which reasonably 
matches the Ci ratio. The Bf ratio at certain lati- 
tudes (no average) are also shown to highlight the 
variation of the effective beams across the sky. 



This 2-column preprint was prepared with the AAS lATgX 
macros v5.2. 
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Fig. 9. — Zoomed version of comparison of maps from Level-S (column # 1) and FEBeCoP (column # 2) at 
two different locations — the ecliptic North pole ( row # 1 & 3) and equator (row # 2 & 4). The difference 
between the Level-S maps and the synfast generated symmetric Gaussian smoothed maps (column # 3) 
show that the effect of asymmetric beams is significant. Finally the difference between Level-S and FEBeCoP 
is very small (column # 4), showing that FEBeCoP simulates the realistic effects as accurate as the detailed 
Level-S simulations. 
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